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Abstract 

An all speed scheme for the Isentropic Euler equation is presented in 
this paper. When the Mach number tends to zero, the compressible Euler 
equation converges to its incompressible counterpart, in which the density 
becomes a constant. Increasing approximation errors and severe stability 
constraints are the main difficulty in the low Mach regime. The key idea 
of our all speed scheme is the special semi-implicit time discretization, in 
which the low Mach number stiff term is divided into two parts, one be- 
ing treated explicitly and the other one implicitly. Moreover, the flux of the 
density equation is also treated implicitly and an elliptic type equation is 
derived to obtain the density. In this way, the correct limit can be captured 
without requesting the mesh size and time step to be smaller than the Mach 
number. Compared with previous semi-implicit methods [11, 13,27], non- 
physical oscillations can be suppressed. We develop this semi-implicit time 
discretization in the framework of a first order local Lax-Friedrich (LLF) 
scheme and numerical tests are displayed to demonstrate its performances. 

AMS subject classification: 65M06,65Z05,76N99,76L05 

Keywords: Low Mach number; isentropic euler equation; compressible flow; 
incompressible limit; asymptotic preserving; Lax-Friedrich scheme. 
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1 Introduction 



Singular limit problems in fluid mechanics have drawn great attentions in the past 
years, like low-Mach number flows, magneto-hydrodynamics at small Mach and 
Alfven numbers and multiple- scale atmospheric flows. As mentioned in [17], the 
singular limit regime induces severe stiffness and stability problems for standard 
computational techniques. In this paper we focus on the simplest Isentropic Euler 
equation and propose a numerical scheme that is uniformly applicable and effi- 
cient for all ranges of Mach numbers. 

The problem under study is the Isentropic Euler equation 

dfPe + V- (p e u e ) =0, 

d f (p e u e )+v(p e u e ®u e ) + ^V/7 e = 0. (1) 

where p e ,p e u e is the density and momentum of the fluid respectively and £ is 
the scaled Mach number. This is one of the most studied nonlinear hyperbolic 
systems. For standard applications, the equation of state takes the form 

p(p)=Ap 7 , (2) 

where A, y are constants depending on the physical problem. 

It is rigorously proved by Klainerman and Majda [15, 16] that when £ — > 0, 
i.e. when the fluid velocity is small compared with the speed of sound [3], the 
solution of £[]) converges to its incompressible counterpart. Formally, this can be 
obtained by inserting the expansion 



p e =P0 + £ 2 p (2 ) + ---, 
U e = U + £ 2 U( 2 ) + • • ■ , 

into (OQ) and equate the same order of £. The limit reads as follows [15, 18]: 



(3) 



P = Po, (4a) 
V-u = 0, (4b) 
<9 / uo + V(uo®u )+V^o = 0. (4c) 

Here po is a scalar pressure that can be viewed as the Lagrange multiplier which 
enforces the incompressibility constraint. Physically, this limit means that in slow 
flows (compared with speed of sound), the factor l/£ 2 in the momentum equation 
in front of the pressure gradient generates fast pressure waves, which makes the 
pressure and therefore, the density, uniform in the domain [23, 24]. 
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For atmosphere-ocean computing or fluid flows in engineering devices, when 
£ is small in standard numerical methods become unacceptably expensive. 
Indeed, £0 has wave speeds of the form 

A =U e ±^vV(p e ), 

where p'(p £ ) is the derivative with respect to p £ . If a standard hyperbolic solver is 
used, the CFL requirement is At = O(eAx). Moreover in order to maintain stabil- 
ity, the numerical dissipation required by the hyperbolic solver is proportional to 
If \X\ = O(^), in order to control the diffusion, we need to have Ax = o(e r ), 
where r is some appropriate constant. Thus the stability and accuracy highly de- 
pend on £. 

Our aim is to design a method whose stability and accuracy is independent of 
£. The idea is to find an asymptotic preserving (AP) method, i.e. a method which 
gives a consistent discretization of the isentropic Euler equations (OQ) when Ax, At 
resolve £, and a consistent discretization of the incompressible limit © when 
£ — > (Ax, At being fixed). The efficiency of AP schemes at the low Mach number 
regime can be proved similarly as in [9]. The key idea of our all speed scheme is 
a specific semi-implicit time discretization, in which the low Mach number stiff 
term is divided into two parts, one part being treated explicitly and the other one 
implicitly. Moreover, the flux of the density equation is also treated implicitly. For 
the space discretization, when £ is 0(1), even if the initial condition is smooth, 
shocks will form due to the nonlinearity of the div(p e u e <S>u e ) term and shock 
capturing methods should be employed here. 

In the literature, lots of efforts have been made to find numerical schemes for 
the compressible equation that can also capture the zero Mach number limit [1,6, 
11,23,24]. In [1], Bijl and Wesseling split the pressure into thermodynamic and 
hydrodynamic pressure terms and solve them separately. Similar to this approach, 
the multiple pressure variable (MPV) method was proposed by Munz et al. in 
[23,24]. There is also some recent work by J. Hauck, J-G. Liu and S. Jin [11]. 
Their approach involves specific splitting of the pressure term. We avoid using 
this splitting, the proper design of which seems very crucial in some cases. 

Some similar ideas can be found in the ICE method, which is designed to 
adapt incompressible flow computation techniques using staggered meshes to the 
simulation compressible flows. The method was first introduced by Harlow and 
Amsdanin 1965 and 1971 [12,13] and is called Implicit Continuous-fluid Eulerian 
(ICE) technique. It is used to simulate single phase fluid dynamic problems with 
all flow speeds. They introduce two parameters in the continuity equation and the 
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momentum equation to combine information from both previous and forward time 
steps. However this method is not conservative, which leads to discrepancies in 
the shock speeds. Additionally it suffers from small wiggles when there are mov- 
ing contact discontinuities. The first problem was solved by an iterative method, 
for example SIMPLE [25], or PISO [14]. In some recent work, Heul and Wessel- 
ing also find a conservative pressure-correction method [27]. All these methods 
are based on the so called MAC staggered mesh in order to be consistent with 
the staggered grid difference method for the incompressible Euler equations [13]. 
Specifically, if we write the simplified ICE technique presented in [2] in conser- 
vative form, we are led to the semi-discrete framework: 



By substituting the gradient of the second equation of © into its first equation and 
using the results of the first equation ©, p e can be updated by solving an elliptic 
equation which does not degenerate when £ — > 0. 

We use a similar idea in our method. However, we do not use the predictor- 
corrector procedure but we rather discretize the problem in a single step. We use 
standard shock capturing schemes which allows to guarantee the conservativity 
and the desired artificial viscosity. We only use implicit evaluations of the mass 
flux and pressure gradient terms to ensure stability and provide an extremely sim- 
ple way to deal with the implicitness. Additionally, we propose a modification 
of the implicit treatment of the pressure equation. Indeed, using a similar idea 
as in [11], we split the pressure into two parts and put ap(p £ ) into the hyper- 
bolic system. This makes the first system no longer be weakly hyperbolic and 
much more stable. The numerical results show the advantage of our method in the 
following sense: 

• The method is in conservative form and can capture the right shock speeds. 

• The non-physical oscillations [10] can be suppressed by choosing the proper 
value of the parameter which determines the fraction of implicitness used in 
the evaluation of the pressure gradient term. The choice of this parameter 
depends on the time and space step and on the specific problem. 




(5) 



(6) 
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In this paper we only use the first order LLF scheme. Higher order space 
and time discretizations will be subject of future work. The main objective of this 
work is to show that the semi-discrete time discretization provides a framework for 
developing AP methods for singular limit problems. Similar ideas can be extended 
to the full Euler equation and more complicated fluid model and have also been 
used in other contexts such as quasineutrality limits [4, 7] and magnetized fluids 
under stong magnetic fields [5]. 

The organization of this paper is as follows. Section 2 exposes the semi- 
implicit scheme and its capability to capture the incompressible limit is proved. 
The detailed one dimensional and two dimensional fully discretized schemes and 
their AP property are presented in section 3 and 4 respectively. In section 5, how 
to choose the ad-hoc parameter is discussed and finally, some numerical tests are 
given in section 6 to discuss the stability and accuracy of our scheme. The ef- 
ficiency at both the compressible and low mach number regime are displayed. 
Finally, we conclude in section 6 with some discussion. 



2 Time Semi-discrete scheme 

Let At be the time step, t n = nAt.n = 0, 1 , • • • and let the V superscript denote the 
approximations at t n . The semi-discrete scheme for the nth time step is 

n n+l _ n n 

Pe Af p£ +V-(p £ u £ )"+ 1 =0, (7) 

(p£U£r+ ^ (p£Ue) " +div( P X^< + a P (p» e )) + l -^-V P {p^) 

= 0, (8) 

where a is an ad-hoc parameter which satisfies a < l/e 2 . The choice of a de- 
pends on the space and time steps and on the fluid speed. When the shock is 
strong, a should be bigger, which means that the system should be more explicit 
to follow the discontinuity more closely. We discuss the choice of a for specific 
equations of state in this paper and test its effect numerically. It depends on the 
required accuracy, the small parameter £ and the shock amplitude in a sometimes 
quite complex way. 

Rewriting the momentum equation ([8]) as 

(p £ u e r +1 = (p £ u e )"-Arv(pX®u2 + ap(p e 1 )) - A;i^^V J P(p e '+ 1 ) 
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and substituting it into the density equation, one gets 

,2 



1 (YP 

Pe + 1 — At 2 ^ AP(Pe + 1 ) = <Kp"X) (9) 

which is an elliptic equation that can be solved relatively easily. Here 

(p» <) = pi - AtV ■ (pX) + A? 2 V • V (pX ® < + MP*)) • (10) 

The Laplace operator in © can be approximated by V(P / (p")Vpg + ) and © 
becomes 

pT x -a; 2 — p^v. (p'(p £ n )Vp £ w+1 ) = *(p»,n»), (ii) 

Though shocks will form for the original system CZ])®, we always add some 
numerical diffusion terms so that 0(p" ,Ug) is smooth. Then so is Pe +1 . When we 
implement this method, Pg +1 can be obtained from © first and u" is then updated 
by the momentum equation ® afterwards. Therefore, apart from the resolution 
of the elliptic equation (PTTI) . the scheme only involves explicit steps. 

We now show that the scheme ©© is asymptotic preserving. We introduce 
the formal expansion 

Ps to = Poc + £ P(\) W + e2 P(2) to + • • • , 

u»=iig(x)+euj 1) (x) + .... ( j 

In the sequel, the V in the index means that the quantity is independent of space. 
When Ax, At are fixed and £ goes to in ©, we formally have AP(pQ +l ) = 0, 
which implies that pQ +1 is independent of space, where pQ +1 is the limit of p£ +1 
when e — > 0. Thus we have 

P °" c+1 A ~ Pgc +V-(p 0c u r+ 1 =0 (13) 

by equating the 0(1) terms in the density equation ©. Integrating (fT3l) over the 
computational domain, one gets 

^ Po^-PoV = -p£ l Jv. (u ) n+1 = -p£ +1 J da *'< +1 - (14) 

As discussed in [11], for wall boundary condition, periodic boundary condition 
and open boundary condition, (fl4l) gives 

Po, +1 =Poc, d5) 



6 



that is po is also independent of time. Thus (1131) also implies 

V-ug +1 =0. (16) 

Then, by using the fact that the curl of the gradient of any scalar field is always 
zero, the curl of the 0(1) terms of the momentum equation ([8]) becomes 

Vx +Vx V(ugOug) =0. (17) 

Thus j 

U ° A ~ U ° + V (u£® ug) + VpJ 2) =0, (18) 

where p n ^ is some scalar field. 

Equations (TT5T) . (fT6l) . (fT8l are the semi-discretization in time of © and thus 
the scheme ©, ([8]) is consistent with the low Mach number limit £ — > of the 
original compressible Euler equation. This statement is exactly saying that the 
scheme is AP. We can see that, in order to obtain the stability and AP properties, 
it is crucial to treat the flux in the density equation © implicitly. 

Letting U = (p e ,p e u e ) r , we can write ©, © abstractly as 

rjn+l _ jjn j 

+ V-F(U n+ 2)+QU n+1 =0, (19) 



At 



where 



f(ti"+^)-( (P^) n+1 ^ o-f ° 7 °\ rim 

Here P is an operator on p e and U n+l / 2 reminds that the flux is partly implicit and 
partly explicit. 

This semi-discretization gives us a framework for developing AP schemes that 
can capture the incompressible limit. Now we are left with the problem of how 
discretizing the space variable. Because shocks can form, considerable literature 
has been devoted to the design of high resolution methods that can capture the 
correct shock speed. Upwind schemes and central schemes are among the most 
widely used Godunov type schemes [19-21]. 

In the present paper, the hyperbolic operator 

At 
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is approximated by an upwind hyperbolic solver and the stiff l/e 1 factor in front 
of the pressure term is treated implicitly. The implicitness of the density flux 
is treated by combining it with the momentum equation. For simplicity, in the 
present work we only consider the first order modified Lax-Friedrich scheme with 
local evaluation of the wave-speed in the current and neighboring cell. 



3 Full time and space discretization: One dimen- 
sional case 

For simplicity, we consider the domain Q. = [0, 1]. Using a uniform spatial mesh 
with Ax = l/M, M being an positive integer, the grid points are defined as 

xj:=jAx, j = 0,l,---,M. 

The flux and Jacobian matrix of (fl9l) become 

W = ( p^Zm ) • - ( _ u|+ ° V(Pe) i ) . CD 

so, the wave speeds are 

A = u e ± Vap'ipe). (22) 
Let Uj be the approximation of U (xj) and let 

A j+ i(0=max(|A J |,|A i+1 |). (23) 

These are the local maximal wave-speeds in the current and neighboring cells. We 
discretize (fT9b in space as follows: 

L + — '-^ + QiU n+l = 0, (24) 

At Ax J 

where F/±i/2(^ n+3 ) is the numerical flux 

F. + ,(^) = l(F+,(^)+F- + ,(t/'^)) (25) 

and 

\/7+l i An n 1l 



F+(U n+ *] 



(p e u £ <g> + ap{p n £j ) +A" M (p £ u £ yj 



n+1 

(p e u e ®u e )'; +1 + ap{p n £j+l ) -A n j+L (p e u e ) n j+1 



and 



Let 



e 2 2Ax 



q = pu, (26) 

and f( 1 \ F ( 2 ) denote the first and second element of F respectively, we can rewrite 
the momentum discretization in (l24l) as follows: 



q^ 1 = q", - AtD*FV (p«, u») - L^^_ (p^) _ ^J*)) . (27) 
Here 

x _ Uj+l/2~Uj-l/2 

Ax • 

By substituting (1271) into the density equation in (1241) . one gets 

P^ 1 - 4 £ ^2 A ' 2 ( P (P^ 2 ) -2 P (pg+ 1 ) + P ( P y_ 1 2 )) = Z>*(p e ",qS), (28) 
where 

DWeMl) = P " £ -AtD^ l \p" £ X £ ) + -I>U)F^\p n £ Xe) (29) 

is a discretization of (j)(p £ ,u' £ ) in (flOl) . We notice that (1281) is a discretization of 
the elliptic equation (fTTI) . We can update q'g +1 through (1271) afterwards. 

To obtain Pg +1 in (|28T) . a nonlinear system of equations needs to be solved. 
One possible way to simplify it is to replace VP(pg +1 ) by Z y (Pg)Vpg +1 , so that 
the following linear system is obtained: 

p *l ' ~ ' ' 4e^ 2 (^V.Hpgfe - Pit 1 ) -Apii-mV-Pl^) 

=ty(p;,<£)- (30) 

This is a five point scheme which is too much diffusive, especially near the shock. 
One possible improvement is that instead of (1301) . we use the following three 
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points discretization 

Pe > ^ (Pe./+l)(P e ;+l-P e y ) ~ P KPejAPej ~ Pej-l) ) 

= D^p'^ £ ). (31) 

After obtaining Pg +1 , we can substitute it into (1271) to get q'^ 1 . 

To summarize, three schemes are proposed here: ( 1281 ), (1271) ; ( 1301) , ( [271 ) and 
(l3"TI) . ((27l) . To investigate the AP property, we take (l30l) . (1271) as an example. 
The proofs for the other two schemes are similar. By substituting the following 
expansion 

Pej = P'oc + £ 2 P( 2)J + ■ ■ ■ > <&j = qoc + eq( 2 ); + ■ • • > (32) 

into (l30l) . the O(^) terms give that p" Q ^j = pQ C +1 is constant in space by using the 
periodic boundary condition, and thus: 

n n+l _ n n+l , 2 n+1 _, 

Pej -P(0)c +£ P(2); + 

Summing (1301) over all the grid points, one gets 

Poc +1 =Poc = Po c , (33a) 

which implies that po is independent of time and space. Thus, the 0(1) terms of 
([301) are 

ApI) {p$ +2 - ) - p'(poc) {p$ - p^U) = o, 

by recalling that the 0(1) terms of both p" and q" are constant in space. Then the 
periodic boundary condition gives 

P(2»=P(& (33b) 
which gives that p"^ 1 is also independent of space. Therefore from ©, (1271) . 

q^ 1 = qoj = q'oc- (33c) 

In one dimension, (1331) is the discretization of (fT5l) . (1T61) . (1T8T) when periodic 
boundary conditions apply and thus is consistent with the incompressible limit. 
In fact all the three methods proposed here are AP. 
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4 Full time and space discretization: Two dimen- 
sional case 



We consider the domain Q. = [0, 1] x [0, 1]. For M\,M 2 two positive integers, we 
use a uniform spatial mesh Ax = 1 /M\ , Ay = 1 /M 2 . The grid points are 

( Xi ,yj) := (iAxJAy), i = 0,- - ,MrJ = 0,--- ,M 2 
Now U = (pE,(^e\(l^) T and f/ !i7 - is the numerical approximation of U(xi, yj). 

Let 



PeU e i 

G X {U)= | p e u 2 + ap(p £ ) | , G 2 (U) 
p e uiu 2 



PeU £ 2 
PeUiU 2 

p £ xx\ + ap(p e 



(34) 



and 







o^p o o 
a^p o o 



Eq. (1191 ) can be written as 

d t U + d x Gi(U) + d y G 2 (U) +QU = 0. 



Denote 



(p e u el )" +1 



G 1 (t/" + 2)= p»(u" £l ) 2 + ap(p n £ 



G 2 (U n+ 2) 



(p£U e2 ) 



n+1 



Pe«el«2 



G 



P>el U £2 

p*(u« 2 ) 2 + a/?(p«; 



\ 

kza<t D xp 

\ lio^ / 



Dfiu 



Uij+l — 



H+lj — Ui-ij 



2Ax l} ~ 2Ay 

Now the eigenvalues of the two one-dimensional hyperbolic equations are 

=m,ui± ^ap'(p £ ), A (2) =u 2 ,u 2 ±v / ap / (Pe)- 
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The fully discrete scheme for the two dimensional problem is 

+^G 2 (t/" +1 / 2 ) + \ (A i j_^D y i j_ -A u+h E? iJ+ )U" + QUf 1 = 0, 



(35) 



where 



and 



l i ~ Ax ' ^J+ U) - Ax 

y _ Ujj — Ujj-l v Ujj+1 — Ujj 

u ij- u - Ay ' Lt ij+ U ~ Ay > 

A^^ma^lA^UA^JJAfUASj}, 
A;; j+r max{|A^|,|A^ 1 |,|Af|,|A^ 1 |}. 



(36) 



Let q e be like in d27l) . Like in one dimension, we can substitute the expressions 
of q"fj ,q2^ into the density equation and get the following discretized elliptic 
equation, 

- (1 ~4p At2 (^( P (Pe&,j) ~2P (pff) +^(P^ 2J )) + (37) 
^2 " 2 J P(P£ 1 ) + P(P^J- 2 ))) = ^(Ps, 

where 

O0y(p e w ,q? e ,q^) 
= p^-Ar(^ 1+ D^ 2 

+A ? 2 (^(p^K 1 ) 2 + ap(p^)+44(p^K 2 ) 2 + ap(p^) 



(38) 
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After obtaining 1 by ( |37l) , q"^\ q'^ 1 can be updated by the momentum equa- 
tion afterwards. 

Similar to the one-dimensional case, the modified diffusion operator using a 
reduced stencil is as follows: 

Peij Ar £ 2 x 

x (^2 { P '(PeiJ + l) (PS/+1 " Pet!) ~ P '(PZ) (P$J - P$J-l) 

+^( P '(Pei + lj)(Pet\j-pZ 1 ) -P'iPZMu-PeAj)) 

= $(Pe^!e,<& £ )- (39) 
Now we prove the AP property of our fully discrete scheme. Here only well- 
prepared initial conditions are considered, which means that there will be no shock 
forming in the solution. Then a can be chosen to be to minimize the introduced 
numerical viscosity. Assuming that the expansions of p e ,u e in (fT2l) hold at time 
t n , when £ — > 0, the 0{\) terms of (O give 



^2 (^(poVi)(pot/ + i -Pot!) -?(fa)(p&-p$j-i\ 

' (^ , (PoVi J 0(Pot\ J -Pot/)- p, (Pou)(Pot 1 -Po i t 1 i J )) =0. 



Ay 



When using periodic boundary conditions, one gets Po ( -t = Po c fr° m ©• The 
time independence of Pg , similar to the one dimensional case, can be obtained 
by summing (1391 over all the grid points. Accordingly we have 

p , ff=PSc + e 2 P'Z)! J + --- (40) 
To prove the limiting behavior of u e i,u £ 2, we do not want to use the density 
equation because the diffusion operator with reduced stencil does not allow us to 
find the corresponding density equation. Therefore, we consider the 0(1) term of 
439]), 

p'^ l -At 2 x 

^ { p '(p n oc) (ps:, +1 - p&y - p, (poc) (ptih - p^-i) 

HpoAA)- (41) 
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Moreover, noting the fact that 



D y ii P(pr 1 )=£ 2 DUp?+ ) l P / (pS c )) +o(s 2 ) 



'(2) 

'\2) 

= E 2 Dfj{pti) Xp, {Pl))+o{E 2 ) 

and similarly, 

the 0(1) terms of the momentum equations of (1351 ) become 

Ar J vp i 'A p / 2 l 'J- <+?j 'J 



Po / J v po 



Comparing (SB with At * (Df- tfESaj) - />: .il42hl) ) , one gets 

D&qg" 1 +Z>;' 7 q'^ 1 = 0(AxAf), (43) 



which is an approximation of (116l) . Moreover, it is obvious that (1421) is a dis- 
cretization of (fTST ). Thus we obtain a full discretization of © in the limit e — > 0. 
Therefore, the two-dimensional scheme is also AP. 



5 The ad-hoc parameter 

In this section we illustrate how to choose At and the parameter a by considering 
the simple state equation P(p £ ) = p £ . In this context, the fully discrete scheme 
(1241) can be written as 

Ei £e _i_ v • a n+l — V • a" 4- V ■ a" — 

qj_-qg _l XT( n n„n^„n j_ „ n ti\ j_ 



At 



+ V(pX®u'e + ape) + 1 ^Vp»+ 1 =0. 
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where V is the centered difference while V stands for the difference of fluxes. The 
latter is defined as follows (in one space-dimension for simplicity): 



{V-F{U n ))j 



f j+ i(u")-f jA (u") 

Ax 



where the flux ^)±i/2 is defined as in (l25l) . By substituting V ■ (q'g +1 — q") from 
the momentum equation of (l44l) into its density equation, one gets 

n"+! — n n 1 — cif 2 

At +^-q n £ -At £2 Ap; 7+1 -ArV-V(pX®Ue + «Pe) = 0- (45) 

The 0(Af) terms behave like a diffusion term which suppresses the oscillations 
at the discontinuity. Assuming that we use a first order explicit LLF scheme, 
the diffusions needed to damp out the oscillations in the mass and momentum 
equations are respectively: 

1 ( \u n E | + ^)AxApg, I (|n» I + ^)AxAq'g . (46) 

Here in (|45T) . besides the 0(At) terms, V • q" also includes some numerical dissi- 
pation. By noting 

p"+ l = pi - AtV ■ q* +1 + 0(AlAx), 
the diffusion for p e now is 

(i (|u e |+Va)Ax+ ^) Ap e " + A? A(pX ® <) (47) 
plus some higher order terms. Moreover, the diffusion for q e is 

1 1 (Y£ 2 

-(\u e \ + y/a) AxAq'g + At ^ Aqg (48) 

and some higher order terms. Comparing (|46|) and (PTTl) . (1481) . in order to suppress 
the oscillations at discontinuities we only need to have 

-(\x^\ + Va)Ax+—^—At>-(\u n e \ + -)Ax, 



that is 



he Ax 

At > 2 _ . (49) 



1 + v / «£ 
15 



Moreover the CFL condition for the explicit part is 



At < a ^ (50) 

max{|u e | + v«| 

where o is the Courant number which is less than 1. We usually choose o to be 
0.5. Then the parameter a should satisfy 

Ax I /— oAx ,, „, , 

-- i <^<^ r -max{|«» W (51) 

according to (|49| ), (l50l) . Then the following constraint on At should hold if we 
want the scheme to be stable and non-oscillatory 

,. _ n Ax ^ oAx 1 
max{u£} + — <—— + _. (52 
2 At At e 

The reason for the occurence of nonphysical oscillations when a = lies in 
the fact that the diffusion is not large enough. In this case, with a simple reduction 
of At, it is likely that the diffusion can no longer suppress the oscillations. This 
is why we need to introduce a to control the oscillations. But, from the analysis, 
no matter the value of a, as long as it is less than 1 /e 2 , the diffusion can never 
be sufficient when At < \eAx. In summary there is no specific way of choosing 
a < 1 /e that can guarantee that the nonphysical oscillations will disappear in 
any case. For well-prepared initial conditions in the low Mach number regime, 
because there is no shock formation in the solution, it is better to choose a as 
small as possible to get better accuracy, but if strong shocks exist in the solution, 
a should be big enough to suppress the oscillations. This is why the choice of a 
depends on the considered problem. 



6 Numerical results 

Three numerical examples will allow us to test the performances of the proposed 
schemes. In fact, three schemes are proposed in section 3 and 4, for example 
in one dimension: the scheme (124)) without linearizing VP(p'^ +l ) is denoted by 
"NL". We need to use Newton iterations to solve the nonlinear system. When 
VP(p" +1 ) is approximated by P'(p'^)VP(p" +l ), the unknowns become a linear 
system. This scheme is represented by "L". "LD" denotes the scheme with the 
narrower stencil (|3TT) . Here we use well-prepared initial conditions of the form © 
and a = 1 for all the test cases. 
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In one dimension, let the computational domain be [a, b] and the mesh size be 
Ax. The grid points are 

xj = a + (j — I) Ax. 

In the following tables, the L 2 norm of the relative error between the reference 
solutions u and the numerical ones U 

c{u) _ \\U-u\\ L 2 _ j i (L J \U J -u(x J )\ 2 y- 

lluLl £(&l«to)l 2 )* 

are displayed. 

Example 1 P(p £ ) and the initial conditions are chosen as 

p e (x,0) = l, p £ (x,0) = l-£ 2 /2 xe [0,0.2] U [0.8,1]; 

p e (x,0) = l+£ 2 , p e (x,0) = l x 6(0.2,0.3]; 

p e (x,0) = l, p e (x,0) = l+£ 2 /2 xe (0.3,0.7] 

p e (jc,0) = l-e 2 , Pe {x,0) = l x 6(0.7,0.8] 

This example consists of several Riemann problems. Shocks and contact discon- 
tinuities are stronger when £ is bigger. We first check the difference of the three 
schemes (|28l . (1301) and (l3"TI) . The CFL condition for the linearized reduced sten- 
cil scheme (l3"TI) is discussed in (ii) and a fixed Courant number independent of 
£ is found numerically. Compared with the first order ICE method using local 
Lax-Friedrich discretization for ©, the improvement of removing nonphysical 
oscillation of our scheme is shown. We investigate the effect of a for different 
values of £ in (iii). In (iv), when a = 1, we numerically test the uniform conver- 
gence order. Finally, the AP property and its advantages are demonstrated in (v) 
by comparing with the fully explicit Lax-Fridrich scheme for the initial Isentropic 
Euler equation (Q~|). 

When £ = 0. 1, the initial density and momentum are displayed in FigureQ]and 
we can see the discontinuities clearly. 

(i) In this example, we choose £ = 0.8,0.3,0.05 corresponding to the com- 
pressible, intermediate and incompressible regimes. The numerical results 
at T = 0.05 of "NL", "L" and "LD" are represented in Figure H Here At is 
chosen to make all these three schemes stable and diminishing At only will 
not improve much the numerical accuracy. The reference solution is calcu- 
lated by an explicit Lax-Friedrich method [19, 20] with Ax = 1/500, A? = 
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Figure 1: Example 1. When £ =0.1, the initial density and momentum are dis- 
played. 

1 /20000. We can see that all these three methods can capture the right shock 
speed. The results of the three schemes are quite close, which implies that 
the linearization idea does simplify the scheme but the "LD" scheme does 
not really introduce less diffusion. When £ is small, though we can no 
longer capture all the details of the waves, the error is of the order Ax which 
is the maximum information one can expect. Numerically, for different 
scales of £, there is not much difference between these three methods. Thus 
in the following one dimensional examples, we only test the performance 
of the "LD" scheme. 

(ii) Because of the explicit treatment of the flux terms in the momentum equa- 
tion, the stability of the 'LD' scheme can be only guaranteed under the 
following CFL condition 

Ax 

At < amin (53) 

i \Ui\ + y/aP'(p e ) 

Here < o < 1 is the Courant number and is set up at initialization. Con- 
sistently with the fact that these three methods are AP, the Courant num- 
ber does not depend on £. Indeed, below, we numerically verify that a is 
independent of £. For £ = 0.8,0.3,0.05, the numerical Courant numbers 
are displayed in Table \T\ and we can see numerically that the biggest al- 
lowed max{w}^: are close to 1 for all £'s. Therefore, a = 0.9 is enough 
to guarantee stability and is numerically shown to be independent of £. By 
contrast, the explicit local Lax-Friedrich scheme for the original Euler equa- 
tion has a stability condition which becomes more and more restrictive as 
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b) ° 




c) 



Figure 2: Example 1 . When T = 0.05, Ax = 1 /200, At = 1 /2000, the density and 
momentum of the "NL", "L" and "LD" schemes for isentropic Euler equation are 
represented respectively by dashed, dash dotted, and dotted lines. The solid line 
is the reference solution calculated by an explicit Lax-Friedrich method [19,20] 
with Ax = 1/500, At = 1/20000. a): £ = 0.8; b): e = 0.3; c): e = 0.05. Left: 
density; Right: momentum. For all e's, these three lines are so close to each other 
that and are not visible in the figure. 
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£ goes to zero. Thus the CFL condition of the standard hyperbolic solver 
At = O(eAx) is considerably improved. 



£ 


max A 


Ax 


stableA? 


Ax 
Af 




0.8 


4.24 


1/100 


1/340 


3.40 


1.25 


0.8 


6.35 


1/200 


1/970 


4.85 


1.31 


0.8 


6.58 


1/400 


1/2420 


6.05 


1.09 


0.8 


6.70 


1/800 


1/5460 


6.82 


0.982 


0.3 


2.64 


1/100 


1/260 


2.60 


1.02 


0.3 


2.70 


1/200 


1/510 


2.55 


1.06 


0.3 


2.76 


1/400 


1/1000 


2.50 


1.10 


0.3 


2.81 


1/800 


1/2050 


2.56 


1.10 


0.05 


2.43 


1/100 


1/260 


2.60 


0.93 


0.05 


2.44 


1/200 


1/490 


2.45 


1.00 


0.05 


2.45 


1/400 


1/960 


2.40 


1.02 


0.05 


2.46 


1/800 


1/1920 


2.40 


1.03 



Table 1: Example 1. The numerical Courant numbers for different £. Here 
max{A} denotes the maximum of max{Aj} defined in (122b until T = 0.1 for all 
time steps. 

(iii) The classical ICE method even in its conservative form introduces some 
nonphysical oscillations, no matter how small the time step is. These oscil- 
lations cannot be diminished by decreasing the time step. Their amplitude 
becomes smaller as the mesh is refined as long as the scheme is stable. 
In this part we show that our method can suppress these oscillations nu- 
merically by choosing a = 1. When T = 0.01, for £ = 0.8,0.3,0.05, the 
numerical results of both our method with a = 1 and ICE calculated by 
Ax = 1 /200, At = 1 /20000 are displayed in Figure [3] The oscillations are 
more important for the ICE method and smooth away when a = 1 . We can 
see that numerical nonphysical oscillations occur in the results of the ICE 
method when £ = 0.8,0.3, but disappear when £ becomes small. This can 
also be seen from (l47b . (l48b . When £ is small the diffusion introduced by 
the implicitness is bigger. These oscillations also go away as time goes on 
due to dissipation. 

(iv) When a = 1, the relative errors of the "LD" scheme for different Ax, At 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 




U U. I \J.C U.J U.*t U.U U.U U./ U.O U.3 I 

c) 

Figure 3: Example 1. When T = 0.01, trie density and momentum for different £ 
are presented. The solid and dashed lines are the numerical results of our scheme 
and ICE with Ax = 1/200^ = 1/20000 respectively, a) £ = 0.8; b) £ = 0.3; 
c)£ = 0.05. 



at time T = 0.1 are shown in Tabled Here Ajc, At do not need to resolve £ 
and the reference solution is obtained by the explicit LLF scheme calculated 
with a very fine mesh Ax = 1 /1280, At = 1/ 128000. We can see that good 
numerical approximations can be obtained without resolving the small £. 
The convergence order is 1/2 when At /Ax is fixed, uniformly with respect 
to £. This convergence order when there are discontinuities is the same as 
the explicit LLF [21]. We can see from Table [2] that refinement in the time 
step does not improve the accuracy much (provided the Courant number 
is appropriately small, like o = 0.7). Take £ = 0.8 as an example. When 
Ax = 1 /320, in order to obtain stability, At should be less than 1/1920. It 
is demonstrated in Table [2] that the error calculated with Ax = 1/320 does 
not decrease much when At is changed from 1 /2880 to 1/12800. Thus as 
long as the scheme is stable, we cannot use a smaller At to obtain a better 
accuracy. This feature is the same as for standard hyperbolic solvers. 



£ 


Ax 


At 


<Pe) 




ratio 






ratio 


0.8 


1/20 


1/180 


9.739* 10- 






1.197 






0.8 


1/40 


1/360 


5.959* 10- 




1.63 


7.484* 10- 


1 


1.16 


0.8 


1/80 


1/720 


3.467* 10 




1.72 


4.180* 10 


1 


1.31 


0.8 


1/160 


1/1440 


1.985* 10- 




1.75 


2.048*10- 


-1 


1.36 


0.8 


1/320 


1/2880 


1.126* 10- 




1.76 


8.477* 10- 


-1 


1.79 


0.8 


1/320 


1/12800 


1.126* 10- 






8.539* 10- 


-1 




0.05 


1/20 


1/70 


4.679* 10- 


-i 




1.355* 10- 


1 




0.05 


1/40 


1/140 


3.305* 10- 


-i 


1.42 


9.574* 10- 


-1 


1.42 


0.05 


1/80 


1/280 


2.353* 10- 


-i 


1.40 


6.758* 10- 


-i 


1.42 


0.05 


1/160 


1/560 


1.655* 10- 


-i 


1.42 


4.430* 10- 


-2 


1.53 


0.05 


1/320 


1/1120 


1.094* 10 


-i 


1.51 


2.538* 10- 


-2 


1.75 


0.05 


1/320 


1/12800 


6.012* 10 


-4 




9.303* 10- 


-i 





Table 2: Example 1. T = 0.1, the L 2 norm of the relative error between the 
reference solution which is calculated with a very fine mesh Ax = 1/1280, A? = 
1/128000 and the numerical results for different £ with different Ax, At are dis- 
played. 

(v) We emphasize the AP property in this final part. For £ = 0.005, the numer- 
ical results at T = 0.01 with unresolved mesh Ax = 1 /20, At = 1/500 and 
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Figure 4: Example 1. By using the "LD" scheme, the density (left) and momen- 
tum (right) for e = 0.005 at T = 0.01 are represented. The circles are the results 
for Ax = 1 /20, At = 1 /500 and the solid line is calculated with Ax = 1 /2000, At = 
1/5000. 



resolved mesh Ax = 1 /2000, At = 1/5000 are displayed in Figure HI while 
the fully explicit Lax-Fridrich scheme is not stable with the same mesh size. 
We do capture the incompressible limit when Ax, At do not resolve £. 



Example 2: In this example we simulate the evolution of two collision acoustic 
waves by the "LD" scheme and test the convergence. We choose a = 1,£ = 
0.1, Ax = 1/100, A? = 1/1000. Here At is chosen to stabilize the scheme and 
decreasing At alone will not improve much the numerical accuracy. Similar to 
Klein's paper [17], P(p E ) and the initial conditions are chosen as 

P(p e ) = pl A , for* e [-1,1] 

p e (x,0) =0.955 + |(l -cos(2ttx)), u e (x,0) = -sign(jc) vT4(l -cos(2rtx)). 

The initial density and momentum are displayed in Figure [5] 

For e = 0.1, the numerical results of the "LD" scheme at different times T are 
shown in Figure [6[ The initial data approximate two acoustic pulses, one right- 
running and one left-running. They collide and their superposition gives rise to 
a maximum in the density. Then the pulses separate again. This procedure is 
demonstrated clearly in Figure [6] 
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Figure 5: Example 2. When £ =0.1, the initial density and momentum are dis- 
played. 



Example 3 In this example, we show numerical results for the two dimensional 
case. Let P(p) — p 2 and the computational domain be (x,y) G [0, 1] x [0, 1]. Be- 
cause no shock will form in this example, we choose a = and the initial condi- 
tion as follows: 

p(x,y,0) = 1 + £ 2 sin 2 (2n(x +y)), 

qi(jc,y,0) = sin(2n(x-y)) + £ 2 sin(2rc(jc+;y)), 

q 2 (x,y,0) = sm(2n(x-y)) +£ 2 cos(2n(x+y)). 

The initial conditions for £ = 0.8 and numerical results at T = 1 calculated with 
Ax—l /20, At = 1/80 are shown in Figure |7J Numerical tests show that a similar 
CFL condition is required as for the one-dimensional case. When £ = 0.05 at 
time T — 1, the numerical results with an unresolved mesh Ax = 1 /20, At = 1/80 
and a resolved mesh Ax = 1/80, At = 1 /320 are displayed in Figure [8j We can 
see that the results using the coarse mesh are much 'smoother' than the one using 
the refined mesh. In this example the amplitude decay due to numerical diffusion 
cannot be ignored. When a coarse mesh is used, the first order method is known 
to have dissipation. This is mainly due to the numerical diffusion term, which 
smoothes out the solution. This phenomenon not only happens when £ is small 
but also when £ is 0(1). We can also see from Figure [8] that when £ is small, 
Z^pie +D y Y>2e is close to 0. 

As a comparison, the numerical solutions of the incompressible limit © with 
and without numerical viscosity are shown in Figure [9j The latter is obtained by 
a difference method based on a staggered grid configuration [13]. This staggered 
difference method is attractive for incompressible flows, since no artificial terms 



24 




Figure 6: Example 2. When £ = 0.1, the density and momentum of the "LD" 
scheme at different times are represented: a) T = 0.01; b) T = 0.02; c) T = 0.04; 



are needed to obtain stability and suppress the oscillations. Because of the stable 
pressure-velocity coupling, solutions with almost no viscosity can be obtained. 
The viscosity introduced here is of the form ^Ax where A is given by d36l) . We 
can see that the amplitude of the wave decays as time evolves even though the 
viscosity is only 0(Ax). In the limit of £ — > 0, (1421) generates a discretization 
of the incompressible limit with O(Ax) numerical diffusion terms. This is why 
the results for Ax = 1 /20, At = 1/80 in Figure [8] are close to those with viscos- 
ity in Figure [9j When the meshes are refined, less diffusion is introduced and 
the solution becomes closer to the solution with no viscosity. The scheme indeed 
catches the incompressible Euler limit and good numerical approximations can be 
obtained without resolving £, which confirms the AP property that is proved in 
section 4. However we need to take care of the numerical diffusion when coarse 
meshes are used. One possible way to improve this is to use less diffusive shock 
capturing schemes at first order or higher order schemes using the MUSCL strat- 
egy for instance [8, 19,20], or to use staggered grid discretizations. 

7 Conclusion 

We propose an all speed scheme for the Isentropic Euler equation. The key idea is 
the semi-implicit time discretization, in which the low Mach number stiff pressure 
term is divided into two parts, one being treated explicitly and the other one im- 
plicitly. Moreover, the flux of the density equation is also treated implicitly. The 
parameter which tunes the explicit-implicit decomposition of the pressure term 
allows to suppress the nonphysical oscillations. The numerical results show that 
the oscillations around shocks of 0{e 2 ) strength can be suppressed by choosing 
a = 1 . The low Mach number limit of the time semi-discrete scheme becomes 
an elliptic equation for the pressure term, so that the density becomes a constant 
when £ — >• 0. In this way, the incompressible property is recovered in the limit 
£ — > 0. Implemented with proper space discretizations, we can propose an AP 
scheme which can capture the incompressible limit without the need for A?, Ax to 
resolve £. 

In this paper we demonstrate the potential of this idea by using the first order 
Lax-Friedrich scheme with local evaluation of the wave speeds. Though this first 
order method is quite dissipative, we can observe that the scheme is stable inde- 
pendently of £ and that the CFL condition is At = O(Ax) uniformly in £. It can 
also capture the right incompressible limit without resolving the mach number. 
Higher order space discretizations like the MUSCL method [8, 19,20] can be built 
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Figure 7: Example 3. When £ = 0.8, the initial density and momentum (left) 
and the numerical result with Ax = 1/20, A? = 1/80 at time T = 1 (right) are 
represented, a) p e ; b) pi e ; c) p 2e . 
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a) 




Figure 8: Example 3. When £ = 0.05, the numerical result with Ax = 1 /20,A/ = 
1/80 (left) and Ax = 1/80, A? = 1 /32O20ight) at time T = 1 are represented re- 
spectively, a): p e ; b): D x pi e +£»'p2e; c): pi e ; d): p 2e 



Figure 9: Example 3. The numerical results of the incompressible Euler limit 
using Ax = 1 /20, At = 1/80 with (left) and without (right) viscosity at time T = 1 
are represented, a): the first element of the velocity uoi; b): the second element 
of the velocity uq2- 
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into this framework. This is the subject of current work. 

This paper provides a framework for the design of a class of all speed schemes. 
Compared with the ICE method [12, 13] and some recent work by Jin, Liu and 
Hauck [11], the idea is simpler and more natural. This framework can also be 
easily extended to the full Euler equation and flows with variable densities and 
temperatures. These extensions and applications [22] will be the subject of future 
work. 
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